Typeset with REVTeX 

Phase Behavior of a Simple Model for Membrane Proteins 



Preprint 



Massimo G. Noro 

Unilever Research Port Sunlight 

Quarry Road East, Bebington Wirral CH63 3JW, U.K. 



o 
o 
o 

(N 

a* 



o 

CO 



o 
o 



> 

(N 

ON 

o 
o 
o 

-I— > 



i 

C 
O 

o 



X 
S3 



Daan Frenkel 

FOM Institute for Atomic and Molecular Physics 

Kruislaan 407, 1098 SJ Amsterdam, The Netherlands 

(February 1, 2008) 



We report a numerical simulation of the phase diagram of 
a simple model for membrane proteins constrained to move 
in a plane. In analogy with the corresponding three dimen- 
sional models, the liquid-gas transition becomes metastable 
as the range of attraction decreases. Spontaneous crystal- 
lization happens much more readily in the two dimensional 
models rather than in their three dimensional counterparts. 



I. INTRODUCTION 

If one picture is worth a thousand words, recent ad- 
vances in X-ray crystallography are providing the equiv- 
alent of a dictionary. Crystallographers are now solving 
the three-dimensional structure of proteins at the rate of 
one or more per day. A bottleneck is the difficulty of 
growing high-quality crystals for X-ray analysis. 

As the success of protein crystallization depends 
strongly on the physical conditions of the initial solu- 
tion [HJ^], much effort has gone into finding the relation 
between solvent conditions and crystallization behavior. 
Rosenbaum et al. 0] analysed the solubility curves of 
a variety of globular proteins and found that they can 
be made to superimpose when expressed in appropriate 
scaled units. What this suggests is that the phase be- 
havior of many globular-protein solutions obeys a law of 
corresponding states. Specifically, they showed that the 
solid-fluid phase boundary of the proteins in solutions 
can be mapped onto the corresponding curve of a simple 
model system (the hard-core Yukawa potential [Q) with 
short-ranged attractions. Such corresponding states be- 
havior suggests that — for a given class of compounds 
— the solubility boundary is only weakly dependent on 
the details of the interaction potential. 

Far fewer crystals have been grown of membrane pro- 
teins than of globular proteins. The reason is simply 
that it is more difficult to crystallize membrane proteins 
than globular proteins. It would, of course, be very in- 
teresting to know if the "generic" features of the phase 
behavior of quasi-two-dimensional proteins are similar to 
those of globular proteins. It is tempting to start such 
an analysis by looking at the corresponding "minimal" 
model for membrane proteins — namely one of circu- 
lar disks with isotropic, short-ranged attractive interac- 
tions. In this paper, we report a numerical study of the 



phase behavior of such model membrane proteins. In our 
study, we vary the range of the attractive interaction be- 
tween the particles. We find that the general topology 
of the phase diagram is indeed similar to that for three- 
dimensional ("globular") proteins. In particular, we find 
that the liquid- vapor transition is preempted by the freez- 
ing curve for particles with sufficiently short-ranged at- 
tractions. However, quantitatively, there are large differ- 
ences. In contrast to the three-dimensional case which 
still exhibits a well-defined meta-stable liquid-vapor co- 
existence curve below the equilibrium freezing curve, we 
find crystallization, rather than fluid-fluid demixing, in 
the corresponding "membrane-protein" model. 



II. MODEL 

We model the effective interaction between membrane 
proteins using an extension of the well-known Lennard- 
Jones potential, 
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where a denotes the hard-core diameter of the particles 
and e the well depth. The width of the attractive well can 
be adjusted by varying the parameter a: the smaller a 
the longer the range of attractions. Figure (Q) plots this 
potential for the values of a used in this paper. Note 
that as a decreases, the range of repulsions increases as 
well, so that the "effective" size of the particle grows. 
It is, however, convenient to compare the simulation re- 
sults for particles that have the same hard-core diameter. 
To estimate the effective hard-core diameter for a given 
value of a, we separate the potential into an attractive 
v a tt and a repulsive v rep contribution in the spirit of the 
Weeks- Chandler- Andersen method ||10| . We then calcu- 
late the equivalent hard-core size of the repulsive part of 
the potential using the Barker-Henderson criterion pi| : 



T e// 



dr 



I _ e Vrep(r)/k B T 



(2) 



In what follows, we use <J e ff as our unit of length. In 
these units, all our model proteins have the same effective 
diameter. 
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FIG. 1. Interaction potential as defined in EqJll correspond- 
ing to a — 50 (solid line), a = 4.5 (dashed line), and a = 0.1 
(dotted line). Note that as a decreases, attractions become 
longer ranged but the range of repulsions increases as well, so 
that the "effective" size of the particle grows. (We have set 
e=l). 



III. COMPUTATIONAL METHODS 

In order to map out the phase diagram of our two- 
dimensional model for membrane proteins, we have used 
a combination of simulation techniques that we discuss 
briefly in this section. 

A. Gibbs-Duhem integration 



This method was proposed by Kofke [|7||8|] , and is based 
on the integration of the Clausius-Clapeyron equation, 
which expresses the slope of the phase boundary in the 
(P,/3) diagram, 



dP 

~d(3 



Ah 

' f3Av 



(3) 



where j3 = l/fc^T, P is the pressure, Ah and Av are 
the differences in enthalpy and volume (per particle) in 
the two coexisting phases, respectively. To compute this 
slope, two simulations are carried out in parallel: one in 
the liquid phase and one in the solid. The two systems are 
held at the same temperature and pressure, but cannot 
interact with each other; during the runs we measure the 
average density (1/v) and enthalpy h per particle and 
thus determine (dP/d(3). Knowing this slope, we can 
then estimate the location of a neighboring point on 
the (P,f3) coexistence curve. The Gibbs-Duhem method 
is straightforward, but there are several ways in which it 
can be implemented, and there are several subtleties that 
require attention. 

The first is the choice of a good starting point. The 
Gibbs-Duhem integration method allows one to trace out 
the (P,P) coexistence curve once one point on this curve 



is known. There are several ways to select this point. 
In one set of calculations, we started the step-wise in- 
tegration at /3 = (i.e. the infinite temperature limit) 
where the phase diagram approaches that of hard disks. 
Here we have used the fluid-solid equilibrium density gap 
reported by Jaster Jy| as the input of our first Gibbs- 
Duhem simulation. However, as the freezing transition 
of hard disks is itself still not completely characterized 
pjfl , this is not necessarily the best option. In fact in- 
dependent free-energy calculations allow us to start the 
Gibbs-Duhem integration at a finite value of (3 where we 
find a strong first-order liquid-solid phase transition. We 
found this second route more reliable. 

The second point is the choice of the form of the equa- 
tion to be integrated. The Gibbs-Duhem method is not 
self-correcting. This means that small numerical errors 
may cause the computed coexistence curve to diverge 
from the true phase-equilibrium curve. To minimize this 
problem, the right-hand side of Eq. |3| must be a smooth 
function of pressure and temperature, so that simple in- 
tegration schemes can be applied with high accuracy. In 
our case we found that at high temperatures a suitable 
slowly varying function was 



d log P _ Ae + PAv 
d\og(3 ~ PAv 



(4) 



where Ae is the difference in the energy per particle be- 
tween the two coexisting phases. In the hard-particle 
limit the PAv term completely overwhelms the energy 
difference, and the slope of the phase boundary plotted 
in the (log P, log j3) diagram approaches — 1. We have 
verified this in our simulations. In the low tempera- 
ture regime the most convenient differential form of the 
Clausius-Clapeyron equation, in agreement with Kofke 
M and Hagen [Q, was found to be: 



d\og(3P 
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(5) 



Equations |I] and || were solved using a second-order 
predictor-corrector algorithm. As such algorithms are 
not self-starting, we initiated the integration by sup- 
plying the values of the integrand and its slope (when 
known), and using a first-order algorithm for the pre- 
diction of the first point; after two integration steps we 
continued with the desired second order procedure. In 
Figure (||) we collect the phase transition points obtained 
from numerical integration of equations ^ and ^| • 

The third point has to do with spontaneous phase tran- 
sitions during a single phase simulation. In two dimen- 
sions there is much less hysteresis in the solid-liquid tran- 
sition than in three dimensions. As a consequence, it may 
happen that, in a constant-pressure simulation of a rela- 
tively small system (in our case: N = 256), a fluid could 
spontaneously transform into a solid, or vice versa. This 
creates a problem for the Gibbs-Duhem simulations that 
involve constant-pressure studies of state points along the 
solid-fluid coexistence line. To prevent such undesirable 
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(and, on a macroscopic scale, irrelevant) fluctuations, we 
imposed a constraint on the degree of crystallinity of the 
system. The degree of crystallinity was measured using 
a global bond-order parameter |l4J . If during a constant- 
pressure Monte Carlo simulation of the liquid (or solid) 
phase, the value of the bond-order parameter of a con- 
figuration was outside the interval typical for the phase 
under consideration, then the configuration was rejected. 
Typically no more than 0.5% of the configurations were 
rejected during a simulation of any state point. 
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(a) (log P, log (3) integration 



log/3 



ature and pressure. The result of such a calculation is 
a (/x,P)-diagram similar to the one shown in Figure (|2|). 
Two phases are in equilibrium at the point where the two 
chemical potential curves cross. Evaluation of the chem- 
ical potential of the fluid branch is straightforward once 
the equation of state of the fluid is known from low den- 
sities up to the density range of interest. In the present 
case, we did this by performing a large number of NPT- 
Monte Carlo simulations at different state points (see 
Figure @j). We then fitted the numerical data for the 
pressure to a convenient fitting function. In the present 
case, we used an ad-hoc generalization of the so-called 
y-expansion that is often used to describe the equation 
of state of hard-body fluids |L5| : 



f3P 



1 — ap 



1 — ap 



1 — ap 



(G) 



where a, b and c are to be determined from the fit. Upon 
integration of the pressure (|6|) between zero density (ideal 
gas limit) and the density of interest, one obtains an ex- 
plicit expression for the chemical potential: 



(3p(p) = In 



pA 2 b/a - c/a 2 + 1 c/2a 2 + bp 
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(b) (log/3P,/3) integration 
FIG. 2. Phase boundary in the two representations used for 
the Gibbs-Duhem integration scheme: (a) in the high temper- 
ature regime (b) in the low temperature regime. 



B. Thermodynamic Integration 

Thermodynamic integration is the method most com- 
monly used to locate the solid-liquid transition. The 
procedure involves the comparison of the chemical po- 
tentials of the liquid and solid phases at equal temper- 



(7) 

where A is the De Broglie wavelength. In the case of 
the solid phase, we fit the calculated (P,p)-isotherms to 
a simple power law of the form ap 2 + bp + c. Integrating 
between a reference density p* and the present density 
yields the chemical potential 

(3p(p) — lap + b (lnp + 1) — (ap* + 6 In// — c/p*) 
+/3r x (p*)+lnA 2 p*-l , 

(8) 

where f* ex (p*) is the excess Helmholtz free energy per 
particle evaluated at the reference density p* . A method 
that is widely used to compute the free energy of a crys- 
talline solid is the so called "Einstein crystal" method 
proposed by Frenkel and Ladd [[16| , which employs ther- 
modynamic integration of the Helmholtz free energy 
along a reversible artificial path between the system of 
interest and an Einstein crystal. The Einstein system 
is used as a reference since there is a simple analytical 
expression for its partition function, which allows a de- 
termination of the absolute value of its free energy. We 
typically performed a series of 10 NVT-simulations at 
the reference density /?*, switching gradually from the 
Einstein crystal to our system of interest, by modifying 
the value of the coupling constant A, where A = 1 corre- 
sponds to the Einstein crystal and A = to the system 
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of interest. The value of the excess free energy takes a 
simple form 1171]: 



pr ex ( P n = - 

N-l 



--In 
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d+1, 1 

lnp+1 — IniV 

H 2N 2N 



ln27r 



(9) 



where d is the dimensionality, N is the size of the system, 
a is the spring constant of the Einstein crystal, and Uq 
the potential energy of the crystal with all the atoms in 
their lattice positions. The difference between the energy 
of the Einstein crystal and that of the system of interest, 
AU = Ueiu — U, enters in the third term in Eq. ^ and the 
integral in this term is evaluated numerically as explained 
in Ref. il. 



fin 
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FIG. 3. Plot of the chemical potential /iasa function of the 
pressure P for the liquid and the solid phase. The transition 
condition is satisfied at the crossing: at a fixed temperature 
the chemical potentials and the pressure are equal. Here we 
neglect the De Broglie wavelength contribution. The calcu- 
lated equilibrium pressure does not change since one subtracts 
the same quantity, namely In A 2 , from the chemical potential 
of both phases. 

Once the transition pressure is known, from the cross- 
ing of the two chemical potential curves of Figure (ft) , one 
simply reads off the coexisting densities from the (P,p) 
diagram. In Figure (Q) we show two typical equations of 
state for the case a = 50 and a = 0.1 where we have 
collected state points using both NPT-simulations and 
/iVT-simulations (i.e. performed in the Grand Canonical 
Ensemble) . At low temperatures it becomes increasingly 
difficult to equilibrate the fluid system, especially if the 
series of simulations is performed as a gradual compres- 
sion. The occasional formation of high density clusters of 
particles generates locally a highly incompressible fluid, 
and a typical compression move is therefore very unlikely 
to succeed. On the other hand, by keeping the chem- 
ical potential constant, the density can be more easily 



increased by adding new particles. The transition cal- 
culated via the thermodynamic integration route is com- 
pletely consistent with our simulation data, and the tran- 
sition is predicted to fall inside the observed "hysteresis" 
loop. 

We have verified the computed (/i,P)-curves of the 
fluid, by calculating independently the chemical poten- 
tial, using the Widom particle- insertion method |18|, in 
an NVT-simulation. 

/3P 
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(a) a = 50 ft = 2 
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(b)a = 0.1 ft = 2 
FIG. 4. Equation of state, or (P,p)-isotherm, calculated 
at the same temperature for two different ranges of inter- 
actions: (a) estremely short (a — 50), (b) extremely long 
(a — 0.1). The simulation data have been collected by: mea- 
suring the average density in NPT-Monte Carlo simulations 
when compressing a liquid (open circles) or when expanding a 
solid (filled circles) , and (open squares) measuring the average 
pressure in pVT-Monte Carlo simulations. 
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C. Gibbs Ensemble Simulations 

A method that focuses specifically on the location of 
the liquid-vapor coexistence curve is the Gibbs-Ensemble 
technique fq] . Here two simulations are carried out in par- 
allel: one in the liquid phase and one in the vapor. The 
two systems are held at the same temperature and are 
allowed to exchange volume and mass, but the total vol- 
ume and total number of particles of the two systems is 
fixed. This ensures that, at equilibrium, the pressure and 
chemical potential of the two systems are the same. As a 
consequence, the conditions for phase coexistence are au- 
tomatically satisfied. Using this technique we calculated 
liquid-gas coexisting densities for the case of long-range 
attractions (a = 0.1), where the liquid-vapor transition 
is stable and the phase diagram shows a critical point as 
well as a triple point. The coexistence data have been 
collected in Table B. Close to the critical point the free 
energy associated with the formation of the liquid- vapor 
interface becomes very small. As a consequence, the free 
energy cost to create an interface in either box becomes 
small, while the formation of such interfaces is entrop- 
ically favorable. For this reason, just below the critical 
point, vapor- liquid coexistence can no longer be observed 
in a Gibbs Ensemble simulation 119] . Therefore the high- 
est temperature at which the coexistence can be observed 
is not a proper estimate of the critical temperature of the 
system; nevertheless it is possible to estimate it by as- 
suming that the temperature dependence of the density 
difference of the two coexisting phases can be fitted to a 
scaling law J2Q]: 



Pli 



Pgas =A(T- T c )~ 



(10) 



where 7 is the critical exponent (for two dimensional sys- 
tems 7 = 0.125), T c is the estimate of the critical tem- 
perature, and A is a constant determined in the fit. Once 
T c is known, it is possible to estimate the critical density 
p c , by using the law of rectilinear diameters pt|: 



Preprint 

been observed in a number of systems E3] and there is, in 
fact, no theoretical reason to rule out first-order melting 
in two dimensions. 
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(a) a = 50 short-range attr. 
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(b) a = 4.5 intermediate-range attr. 



P a lff 



Pliq ~r Pgas 



= p c + B(T- T c ) 



(11) 



where B is an adjustable parameter. 



IV. RESULTS 

We have mapped the phase boundaries for three dif- 
ferent ranges of attractions corresponding to a = 50, 4.5 
and 0.1. Our results are summarized in Figure (ph. The 
first point to note is that we find clear evidence of a 
first order transition between the solid and fluid phase 
at finite temperatures. This finding would be trivial in 
three dimensions but not in two. In fact there is consid- 
erable evidence that melting in two dimensions may be a 
continuous phase transition pl|]22| ]. On the other hand 
evidence for first order two-dimensional melting has also 
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FIG. 5. Phase diagrams in the (/3, p) representation, shown 
for increasing range of attractions: (a) a = 50, (b) a — 4.5, 
and (c) a = 0.1. The circles refer to the Gibbs-Duhem in- 
tegration method. Here we show (open circles) a series for 
increasing j3 (decreasing temperature) and (filled circles) a 
series for decreasing /3. The squares (open squares) are the 
result of thermodynamic integration, and the line just a guide 
to the eye. The fluid-fluid equilibrium densities are calcu- 
lated with Gibbs Ensemble simulations (open triangles) when 
possible, otherwise we estimated the spinodals (dashed line) 
extrapolating the isothermal compressibility. See text for de- 
tails. 



A. Short-range 

In the first panel of Figure (0) we show the calculated 
phase diagram for the shortest range considered. A rel- 
atively small system size of N — 256 particles was used 
in the mapping of the phase boundary using two par- 
allel NPT-simulations combined into the Gibbs-Duhem 
method. In order to shorten simulation times we have 
truncated, but not shifted, the potential at r = 3.5a, 
and maintained a neighbor list of particles within a ra- 
dius of r = 5. Oct. In order to prevent the solid system 
from melting and the liquid from crystallizing, we have 
used the artificial constraint on a crystallinity order pa- 
rameter (see Paragraph [II). However, we also ran a 
few simulations without this artificial constraint and we 
verified that the phase coexistence was not an artifact 
due to the constraint. We performed two sets of integra- 
tions, one where we increased (3 step-wise (empty circles 
in the figure), and one in the opposite sense (filled cir- 
cles). An independent evaluation of the crystallization 
boundary was obtained with thermodynamic integration. 
The squares in the figure represent the calculated coex- 
isting densities, and are connected with a simple fitting 
function which is only meant to be a guide to the eye. 
Note that the curves for increasing and decreasing j3 do 
not superimpose everywhere. This is due to the limited 
numerical accuracy of our Gibbs-Duhem integration. 

It is interesting to compare our results with the phase 
diagram calculated for particles interacting through the 
same potential and with same range , but in three di- 
mensions. Constraining the system to two dimensions 
causes the "shoulder" in the crystallization curve to be- 
come flatter and to move to lower temperatures. The 
latter effect is not surprising as there are more neigh- 
bors in d = 3 than in d = 2. For instance, a sphere in 
a d = 2 close-packed structure has 6 neighbors and in 
d = 3, 12 neighbors. All other things being equal, in 
three dimensions the freezing temperature is raised by a 
factor proportional to the number of neighbors. Next, 
the effect of thermal fluctuations increases as the dimen- 
sionality of the system is reduced. In d — 1, solids can 
not exist because of thermal fluctuations. In d = 2, there 
is no true long-range positional order, although there is 
a phase transition separating the liquid and the crys- 



talline phase. We should expect that, due to the stronger 
fluctuations in two dimensional fluids, the (meta-stable) 
liquid-vapor critical temperature, T c , should be reduced 
compared to the corresponding three-dimensional model. 
We have attempted to locate the metastable fluid- 
fluid equilibrium by performing Gibbs-Ensemble simu- 
lations in the region where we estimated demixing of the 
metastable fluid phase to occur. However, in these sim- 
ulations we saw no fluid-fluid demixing. Rather, we no- 
ticed a strong tendency towards spontaneous crystalliza- 
tion. This behavior is in contrast with analogous cal- 
culations in the corresponding three-dimensional system 
Q. We argue that in two dimensions the local hexago- 
nal structure of the dense fluid, is similar to that of the 
solid. Only small density fluctuations are necessary to 
overcome a (presumably) small free energy barrier for the 
formation of the critical nucleus. Nevertheless we can still 
provide an estimate of the fluid-fluid metastable equilib- 
rium by extrapolating to the temperature where the in- 
verse isothermal compressibility, k^ 1 = — V (dP/dV) T , 
of the liquid phase vanishes. This is straightforward, 
since the equation of state of the liquid is known for 
several temperatures from the thermodynamic integra- 
tion procedures. See for example Figure (0). The inverse 
isothermal compressibility can be rewritten as: 



/3 



-i fd(3P\ 



(12) 



We assumed that the inverse compressibility depends lin- 
early on temperature (this is not true close to the criti- 
cal point, but there we have no data anyway). The set of 
points where the extrapolated /3k _1 vanishes, provides us 
with our estimate of the meta-stable liquid-vapor spin- 
odal. This estimate is shown as a dotted curve in Figure 

(!)■ 



B. Intermediate- range 

We used thermodynamic integration to map out the 
solid-liquid equilibrium for longer range of attractions. 
Even though the calculation is limited to a few selected 
temperatures, this method has the advantage that it is 
more robust than the Gibbs-Duhem. The curve con- 
necting the squares in Figure (||) is only a guide to the 
eye. As we increase the range of attractions, we expect 
the metastable critical point to move to higher temper- 
atures. Indeed our estimate of the critical point, ob- 
tained using the extrapolation method described above, 
predicts that, for the intermediate range (a — 4.5), the 
spinodal just about touches the crystallization boundary. 
We attempted to perform Gibbs ensemble simulations to 
study the liquid-vapor transition in this model system, 
but again we encountered a strong tendency of the dense 
phase to crystallize. 
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C. Long-range 

The situation changes quite dramatically when the 
range of attraction is increased even more. The liquid- 
gas transition becomes stable and Gibbs-Ensemble simu- 
lations are successful in locating the phase boundaries. In 
Table | we collect the coexistence densities as a function 
of the inverse temperature. The uncertainties Ap quoted 
in the table are not the errors in the average densities, 
rather they refer to the half- width of the histogram indi- 
cating the probability of finding a certain density during 
the simulation. By fitting the liquid and gas densities 
to the law of rectilinear diameters, we extrapolate the 
critical point at T c = 0.418 and p c = 0.134. The fitting 
curve is portrayed in the third panel of Figure (g). The 
open circles are the Gibbs-Ensemble simulation coexis- 
tence densities for the liquid-gas equilibrium, while the 
open squares refer to the solid-fluid equilibrium and were 
determined through thermodynamic integration. At very 
low temperatures, though, calculating the fluid branch 
of the equation of state becomes quite a difficult task, 
because equilibration times become increasingly longer. 
Here we have estimated the crystallization boundary by 
Grand-Canonical (/iVT) simulations of very low density 
systems. 




0.0 2.0 4.0 6.0 8.0 10.0 r ' Ue ^ 

FIG. 6. Pair correlation function g(r) calculated in the 
short-range attraction a — 50 system at T* = 2.5 for different 
densities: (dotted line) typical gas density p* = 0.13, (dashed 
line) liquid system exactly at equilibrium p* = 0.5680, 
and (solid line) the corresponding coexisting solid system 
p* = 0.8267. 



D. Structural properties 

In order to investigate the structural properties of the 
phases identified in the previous sections, we have per- 
formed simulations on a larger system of 450 particles, 
at a single temperature, for the case of a short-range po- 
tential. The use of this larger system enables the calcu- 
lation of the pair distribution functions for large separa- 
tion: this is especially important in order to study the de- 



cay of the bond-order parameter close to the solid-liquid 
transition. Monte Carlo simulations were performed in 
the NPT-ensemble. The structures of the two coexisting 
phases were characterized using the radial distribution 
function g(r) (see Fig. (H)) and the bond-order correla- 
tion function g§ (r) , which correlates over the value of the 
local bond-orientational order parameter ipQ between two 
particles. This is defined by: 



9e(r) 



9(r) 



(13) 



where the local bond orientational order for particle i at 
a position r$ is given by 



tpei^i) 



J2k w(r ik ) e-x.p(6i9 ik )) 

J2k w ( r ik) 



(14) 



In this expression the summation is over the neighboring 
particles k of particle i and Oik is the angle between the 
vector (i"i — r-fc) and a fixed reference axis. We used a 
weighting function w to define nearest neighbors f24|| . In 
the present study, we chose w(r) such that it is unity for 
a separation of r < 1.6(7 and zero for r,fc above 1.8a with 
a linear interpolation between the two endpoints. The 
upper limit of the weighting function was chosen such 
that all the particles included in the first peak of the pair 
distribution function contribute to "06 • The bond-order 
correlation function is large if local bond-order parame- 
ters are correlated over large distances (as they are in the 
crystalline solid). In the isotropic liquid, bond-order cor- 
relations decay exponentially. However, in the hexatic 
phase, an algebraic decay is expected. Our results for 
ge(r) are shown in Figure (Q): it is apparent that g&{r) 
rapidly tends to zero in the liquid phase, and therefore 
the liquid phase at coexistence is not even close to be- 
coming hexatic. Of course, in the coexisting solid phase, 
the bond orientational correlation function tends to a 
non-zero value, as it should. These findings support our 
"thermodynamic" observation that the fluid-solid transi- 
tion in these model systems appears to be first order. 
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FIG. 7. The bond orientational correlation function g&(r) 
calculated for the systems described in Fig. (pi). Note how 
the low density system carries no bond correlations even at 
short distances; the correlation in the liquid is lost within 4 — 5 
particle diameters, but in the solid the bond-order parameter 
tends to a non-zero value in the limit of long particle-particle 
separation. 



V. CONCLUSIONS 

We have studied the phase behavior of a simple model 
system that is meant to mimic membrane proteins con- 
fined in a quasi-two-dimensional geometry. To study the 
phase behavior, we used a variety of complementary sim- 
ulation techniques. Where possible we have used various 
routes to estimate the phase boundaries. Our study fo- 
cused on the influence of the range of the attractive inter- 
actions on the topology of the phase diagram. For long- 
ranged attractions the phase diagram displays a stable 
liquid-vapor critical point and a solid-liquid-vapor triple 
point. As the range of attraction is decreased, the stable 
liquid-gas transition becomes metastable and the critical 
point moves into the solid-fluid two-phase region. Quali- 
tatively, the phase diagrams of the two-dimensional sys- 
tems that we studied are similar to those of their three- 
dimensional counterparts. However, quantitatively, there 
are large differences. Most importantly, we find that 
in systems where the liquid-vapor coexistence curve has 
moved below the freezing curve, there is virtually no bar- 
rier to crystallization. This suggests that membrane pro- 
teins with effectively isotropic interactions should easily 
form two-dimensional crystals. Two questions arise: 1.) 
how do these two dimensional crystals proceed to form 
three-dimensional crystals and 2.) to what extent is the 
ease of 2d crystallization changed by anisotropy in the 
protein-protein interactions? 
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Pgas 



Ap 



Pliq 



Ap 



2.60 


0.015 


0.010 


0.250 


0.005 


2.55 


0.020 


0.010 


0.245 


0.005 


2.525 


0.020 


0.010 


0.240 


0.010 
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2.50 


0.0225 


0.010 


0.235 


0.010 


2.475 


0.020 


0.015 


0.230 


0.015 


2.45 


0.032 


0.015 


0.228 


0.010 


2.425 


0.040 


0.020 


0.230 


0.020 


2.40 


0.060 


0.015 


0.220 


0.015 


2.375 


0.070 


0.020 


0.200 


0.020 


2.35 


0.090 


0.020 


0.180 


0.020 



TABLE I. Coexistence data for the liquid-gas equilibrium. 
/3 is the inverse temperature and p is the number density. 
The uncertainties Ap quoted here refer to the half-width of 
the histogram indicating the probability of finding a certain 
density during the simulation. By fitting p gas and pu q to the 
law of rectilinear diameters, we extrapolate the critical point 
at (3 C = 2.392 and p c = 0.134. 



